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PREFACE 

This  Memorandum  presents  a  method  of  estimating  the  bearing  angle 
of  an  incoming  plane  wave  using  an  arbitrary  ground  array  of  sensors. 
It  was  prepared  for  the  Advanced  Research  Projects  Agency's  VELA 
Analysis  study.  The  project  is  a  broad  and  continuing  system-oriented 
study  of  the  detection  of  nuclear  bursts  above  the  ground. 

The  Memorandum  should  be  useful  to  those  concerned  with  acoustics 
and  seismology,  as  well  as  those  interested  in  data  processing. 
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SUMMARY 

This  study  is  concerned  with  developing  data  processing 
techniques  to  obtain  bearing  angle  estimates  of  plane  sonic  waves 
using  arbitrary  ground  arrays  of  microphones.  The  evaluation  of 
the  accuracy  obtainable  as  measured  by  the  rms  bearing  angle  error 
is  computed  in  detail  for  a  16-station  square  array.  A  novel 
feature  of  the  method  is  that  the  ground  trace  velocity  of  sound  need 
not  be  known  a  priori  or  measured  independently,  but  can  be  derived 
from  the  same  measurements  as  the  bearing  angle. 
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I.  INTRODUCTION 


STATEMENT  OF  THE  PROBLEM 

Given  an  array  of  nondirectional  microphones  which  measure  sound 
pressure,  it  is  desired  to  measure  the  bearing  angle  of  an  arriving 
plane  acoustic  wave  in  the  infrasonic  regiom  ,  that  is,  in  the  frequency 
range  of  from  .  1  to  1  cps.  The  array  may  be  of  arbitrary  geometry  in 
the  ground  plane.  A  novel  aspect  of  the  problem  is  that  the  local 
velocity  of  sound  propagation  is  not  presumed  known  except  for  a 
nominal  value  of  c^  =  344  m/sec.  The  actual  velocity  may  deviate  by 
5  to  10  percent.  The  local  ground  trace  of  sound  propagation  is  also 
obtainable  from  the  measurements  as  described;  however,  estimates  of 
the  elevation  angle  of  the  plane  wave  are  not. 

DESCRIPTION  OF  THE  MODEL 

The  plane-acoustic  wave  is  presumed  to  be  generated  a  large 
distance  from  the  array.  As  the  sound  wave  is  propagated  through  the 
atmosphere,  the  wave  undergoes  changes  in  both  orientation  of  the  phase 
plane  and  amplitude.  The  amplitude  decreases  slightly  due  to  atmos¬ 
pheric  absorption,  but  primarily  due  to  the  dilution  of  the  sound 
energy  over  a  greater  volume.  Superimposed  on  these  systematic  effects 
there  are  also  random  changes  in  phase  at  each  point  on  the  phase  plane 
caused  by  turbulence  in  the  atmosphere.  Thus,  the  wave  which  arrives 
at  the  array  is  not  strictly  a  plane  wave.  The  surfaces  of  constant 
phase  are  taken  to  consist  of  a  plane  plus  random  deviations  from  the 
plane.  An  excellent  discussion  of  the  propagation  properties  of 
infrasonic  sound  waves  through  the  atmosphere  is  given  in  Ref.  1. 


The  acoustic  plane  wave  energy  (noted  as  the  signal)  is  presumed 
to  be  small  compared  to  the  atmospheric  turbulence  pressure  (noted  as 
noise)  in  the  same  frequency  range.  It  is  assumed  that  the  signal  has 
been  detected  by  other  means  and  that  the  gross  direction  (within, 
say,  one  quadrant)  of  the  wave  has  been  determined.  This  Memorandum 
is  therefore  not  concerned  with  the  detection  problem  but  with  the 
improvement  of  the  estimate  of  the  local  bearing  angle  of  the  sound 
wave.  The  data  at  each  array  point  are  the  result  of  processing  the 
received  data  through  noise- reducing  line  microphones  to  improve  the 
signal- to-noise  ratio.  This  Memorandum  does  not,  however,  attempt  to 
evaluate  the  nature  of  the  background  noise  or  the  effects  of  various 
data  processing  operations  on  the  statistical  properties  of  the  signal 
and  noise.  These  problems  will  be  considered  in  future  studies.  A 
class  of  bearing  estimation  methods  are  developed  and  the  effect  of 
two  spec  i  i  ic  mettiods  is  evaluated  for  certain  standardized  error  models. 
The  measure  of  merit  used  is  a  normalized  standard  deviation  of  bearing 
angle  error,  noted  as  L(0)  .  A  set  of  computations  of  L(0)  is  performed 
for  a  square  array  consisting  of  16  equally  spaced  array  points. 


II.  DISCISSION  OF  THE  ESTIMATION  METHOD 


The  basic  data  required  are  the  transit  time  of  the  wave  from  a 

fixed  station  or  array  point  with  coordinates  (x  ,  y  )  to  each  of  the 

o  o 

other  stations  with  coordinates  (x.,  y,).  Let  this  Lime  be  noted  as 

i  i 

T..  Then  the  estimation  process  involves  the  following:  If  each  of 
the  values  of  t.  is  plotted  in  the  (X,  Y)  plane 


where  d  is  a  normaliEing  jCale  factor,  at  the  value  (X^ ,  Y^),  it  will 
be  shown  that  the  transit  time  for  a  plane  wave  can  be  represented  as 

-  =  A,X  +  AY  +  A~XY  +  A.X2  +  AJ2 
1  2  3  4  5 

The  bearing  angle  9  and  the  ground  velocity  c  are  estimated  from  the 
coefficients  A^  and  The  process  then  involves  estimating  the 

coefficients  A^  by  curve  fitting  of  Eq .  (1)  to  the  data  set 
I*  },  i  *  1,  2,  .  ..,  N  for  an  (N  +  1)  station  array.  Let  the 
measured  value  ol  ^  be  given  by 

T .  =  ^  .  +  AT  . 

it  t 

where  t  is  the  "true"  transit  time  given  by  Eq .  (1)  and  At.  is  the 

t  1 

random  error  in  transit  time  due  to  such  causes  as  initial  phase 
errors  or  deviation  from  the  plane  phase  surface,  and  errors  in 
estimating  t.  from  the  processing  of  signals  from  the  array  micro¬ 
phones.  As  an  example,  an  obvious  method  of  estimating  is  by 
cross  correlation.  That  is 
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+t 

T .  =  p  (t)  =  lim  77  I  s  (u)  s . (u  +  T )  du  (3) 

l  max  t  2t  J  o  i 

t  —  00  - 1 

t  h 

where  s.  is  the  observed  signal  from  the  i  array  point  and  s^  is 

the  observed  signal  from  the  reference  array  point.  The  observed 

and  the  statistics  of  At.  will  be  determined  by  such  factors  as 

i 

the  actual  interval  of  time  over  which  the  cross  correlation  is  performed 
whether  the  computation  of  the  cross  correlation  is  for  sampled  data 
or  for  continuous  data,  and  the  sampling  rates;  the  time  and  space 
correlation  properties  of  both  the  signal  and  noise  components  of  the 
observed  signal;  and  the  difference  in  the  initial  timing  errors  of 
each  observed  signal  due  to  initial  phase  errors  of  the  plane  wave. 

The  above  effects,  as  well  as  alternate  methods  for  generating 
the  t ^ ,  will  be  considered  in  subsequent  studies.  For  the  purpose  of 
this  study,  the  random  variable  At^  is  assumed  to  have  the  following 
propert ies 


E(AT.)  =  0 

2 

E(At.  At .)  =  c  6ii  Case  A 
^  J  J 

2 

E(At.  ATj)  =  a  Case  B 


(4) 


where  E(  )  signifies  the  expec 
i  /  j .  The  quantity  ijf  is  a 
is  assumed  to  have  the  form 


ted  value  and  6^  =  1 ,  i  =  j ;  =  0 , 
normalized  correlation  coefficient 


and 


♦ij  =  cxp  {-rij/k} 


(5) 


,  ,  .  th  ,  .  th  . 

where  r  is  the  normalized  distance  between  the  l  and  j  station 

ij 


1 


5 

/  2  r 

r.  .  =v(x.  -  X.)  +  (Y.  -  Y.) 
ij  i  J  i  J 

and  k  is  a  constant  s  0. 

The  coefficients  A  of  Eq .  (2)  are  obtained  by  generalized 
least-squares  procedures  as  follows:  Let  N  +  1  be  the  number  of 
stations  so  that  the  number  of  transit  times  t.  measured  from  the 

l 

reference  station  is  N.  Define  an  N  X  N  positive,  definite,  symmetric 
matrix  p  with  elements  p  }.  Then  let 

r  ‘•'UjV 

N  N 

o  *  1  i  <Ti  -  V  Pij<Tj  -  Tj>  <6> 

i=l  j-1 

The  values  of  A.  are  selected  which  minimize  Q. 

J 

When  the  following  conditions  hold,  the  solutions  are  as  indicated: 


p  =  I  (identity  matrix)  - 

least-squares  solution 

(7a) 

p  =  diag  {p  }  =  weighted 

least-squares  solution 

(7b) 

p  =  fp  }  =  generalized 

1  u,vJ 

weighted  least-squares  solution 

(7c) 

pijf  =  I  =  minimum  variance 

solution 

(  7d) 

The  general  formulation  for  Eq.  (7c)  is  shown,  from  which  Eqs.  (7a), 

(7b),  and  (7d)  are  given  as  special  cases.  Computations  of  L(0)  are 
performed  for  the  square  array  consisting  of  16  equally  spaced  arrays 
separated  by  distance  d  between  x  and  y  coordinates  of  adjacent  stations. 
Similar  computations  are  performed  for  the  linear  case  (A^  =  A^  =  A^  = 

0)  and  for  certain  subsets  of  stations  to  measure  the  improvement 
rate  in  L(Q)  as  more  stations  are  processed. 

CORRELATION  MATCH  VERSUS  MISMATCH 

The  effects  of  mismatching  the  weighting  matrix  p  and  the  At^ 


* 
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correlation  matrix  '•  are  computed  as  follows 

Case  I:  p  =  I;  ,  given  by  Eq.  (4),  Case  B;  k  >  0 

Case  II:  p  -  ,-l;  ,  given  by  Eq.  (4),  Case  B;  k  >  0 

That  is,  the  data  correlation  is  actually  as  given  by  Eq.  (5), 

but  a  least-squares  solution  Eq.  (7a)  is  used.  Note  that  as  k  -  0 , 
min(ri.)  fixed,  i  *  j ,  ■  .  I,  so  that  the  solution  for  the  A.  approaches 

the  matched  condition  given  by  Eq.  (7d),  i.e.,  Case  II.  The  matched 

A 

condition  is  optimum  in  the  following  sense.  The  estimates  A.  obtained 

are  random  variables  with  zero  mean  and  covariance  matrix  B  =  fB  ]  u.v 

L  u,vJ 

1,2, ...5.  B  is  positive  definite  (in  the  quadratic  case,  A^ ,  A,, 

A^  1  0)  and  =  the  covariance  matrix  of  any  other  linear  unbiased 
estimator  of  A '  =  (A^ ,  ,  ...  A,.). 

Thus,  a  comparison  of  the  values  of  L(0)  for  the  matched  and 
mismatched  case  shows  how  much  is  gained  by  using  a  minimum  variance 
estimator  as  opposed  to  a  least-squares  estimator.  Comparison  of  the 
subsets  N  =  3,  7,  15  (linear)  and  N  =  7,  15  (quadratic)  shows  how 
much  is  gained  by  using  the  additional  stations.  Finally,  a  comparison 
of  L(9)  for  the  quadratic  curve  fit  and  the  linear  curve  fit  shows 
how  much  additional  root  mean  square  error  is  caused  in  assuring  an 
unbiased  estimate  of  the  bearing  angle  0.  It  may  be  desirable  to 
accept  a  linear  model  for  Eq.  (1)  and  a  small  bias  in  e  with  smaller 
rms . 

ADVANTAGES  OF  GENERALITY  OF  METHOD 

The  technique  does  not  depend  on  the  specific  geometry  of  the 
array.  Thus,  the  method  lends  itself  to  field  data  measurement 
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procedures  since  dropping  bad  data  does  not  upset  Lite  computations. 
Further  unreliable  data  can  bo  weighted  to  have  less  effect. 

In  Appendix  A  the  solution  for  the  coefficients  A.  is  given  in 
terms  of  the  observations  T. .  The  equations  for  the  covariance  matiix 
II  required  to  evaluate  the  variances  of  0  and  c  are  also  derived. 

In  Appendix  B  he  justification  for  Eq.  (1)  and  the  interpretation 
of  the  coefficients  in  terms  of  the  geometry  of  the  plane  wave  and 
local  meteorological  conditions  are  shown.  The  condition  for  accepting 
a  linear  model  is  derived;  that  is,  setting  A^  -  A^  =  A,.  =  0  in 

Eq.  (1). 

In  Appendix  C  the  normalized  rms  bearing  angle  error  L(0)  is 
derived  in  terns  of  the  covariance  matrix  B  of  the  parameter  estimates 

A 

A.  The  ground  trace  transit  time  to  travel  the  distance  d  given  by 
the  scale  factor  in  defining  X  and  Y  is  defined  as  tq.  The  normalized 
rms  error  in  T  ,  M(Q) ,  is  also  derived  in  terms  of  the  same  variables. 
The  results  are  presented  in  tables  following  Appendix  D.  Table  1 
presents  L(e)  for  the  Case  I,  (p  =  I)  versus  selected  values  of  k 
for  the  linear  case  0  =  0°,  15°,  30°  and  45  and  N  -  3,  7,  15.  lhe 
value  of  the  ratio 

R-  L<e)u.s ./ 

is  also  shown  in  the  table  where  L(c.)  .  is  the  matched  processing 

min 

case  p  ^  =  I.  The  value  of  R,  which  is  =  1,  shows  the  gain  obtained 

by  using  the  matched  processing.  The  same  information  is  presented 

for  the  quadratic  case  for  N  =  /,  15  in  Table  2.  In  Table  3,  the 

o 

same  information  is  presented  for  M(0)  for  0  =  45  . 


As  shown  in 
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Appendix  C,  M(0)  =  L(Q)  for  0=0°  and  90°  and  the  maximum 
|M^(9)  -  L^(0)  |  occurs  at  9  =  45°. 

In  Appendix  D  the  computations  of  M(0)  and  L(0)  for  a  specific 
square  array  of  sensors  is  described.  The  Fortran  program  for  the 
computations  is  given.  The  results  are  presented  in  figures  following 

Appendix  D.  Figures  3  to  8  are  plots  of  L(0)  versus  0  for  the  para¬ 

meters  as  plotted  for  fixed  values  of  k  =  .125,  4  and  256.  In  Fig.  3, 
k  =  .125  is  taken  as  indicating  independent  timing  errors  so  that, 
since  =  I ,  the  least-squares  solution  is  a  matched  solution.  For 
Fig.  4,  k  =  4  is  taken  as  a  moderately  mismatched  leas t- squares 
solution.  In  Fig.  5,  k  =  256  is  taken  as  a  heavily  mismatched  solution. 
For  Figs.  6,  7  and  8,  the  matched  solution  is  presented  for  the 
corresponding  cases  of  k  of  Figs.  3,  4  and  5.  Figures  9  and  10  show 
L(@)  versus  k  for  fixed  0,  for  the  linear  case  N  =  3,  7  and  15  and 
the  quadratic  case  N  =  7,  15.  Figure  9  is  for  0=0°  and  Fig.  10  is 
for  0  =  45  .  Both  are  for  Case  I,  p  =  I.  The  same  data  are  presented 

in  Figs.  11  and  12  for  Case  II,  the  matched  case  for  p  ^  =  I.  Other 

angles  are  obtainable  from  Tables  1  and  2.  Note  that  for  Case  I, 

L(0)  is  labeled  L(0)  and  for  Case  II,  L(0)  .  .  The  value  R  of 

L.i>.  mm 

Tables  1  and  2  is  given  by 


L(0) 


R  =  - 


L(0)  . 

min 


L.  s .  > 


where  corresponding  values  of  each  of  the  parameters  are  used  in  the 


ratio . 


By  inspection  of  the  tables  and  graphs  conclusions  can  be  made 
as  to  the  accuracy  in  bearing  angle  obtainable  as  a  function  of  bear¬ 
ing  angle  0,  increasing  station  numbers,  using  linear  versus  quadratic 
curve  fitting,  the  degree  of  mismatch  for  the  least-squares  estimate, 
and  the  accuracy  gain  using  a  minimum  variance  estimate. 

For  example,  in  Fig.  4  for  linear  curve  fitting  there  is  apparently 
little  to  be  gained  at  any  angle  0  by  processing  more  than  N  =  3. 
However,  in  the  quadratic  case  there  is  a  substantial  gain  by  going 
from  N  =  7  to  N  =  15.  This  gain  is  dependent  on  0  and  increases 
monotonically  from  0=0  to  0  =  45  . 
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III.  CONCLUSIONS 

A  method  of  estimating  the  bearing  angle  of  a  plane  sonic  wave 
using  an  arbitrary  ground  array  of  sensors  has  been  developed.  The 
method  does  not  require  knowledge  of  the  propagation  velocity  of 
sound.  In  fact,  the  ground  trace  velocity  of  sound  can  be  derived 
from  the  data  processing. 

Equations  for  evaluating  the  rms  bearing  angle  error  and  the  rms 
ground  trace  timing  error  were  developed. 

Computations  of  L ( 0)  and  M(9) ,  the  normalized  rms  errors,  were 
performed  for  a  specific  square  array  consisting  of  16  equally  spaced 
microphones.  For  this  array,  the  computations  demonstrate  the  accuracy 
obtainable  in  terms  of  the  rms  timing  errors  and  provide  a  basis  for 
determining  how  to  efficiently  process  the  field  data. 

«v 

Subject  to  certain  mild  restrictions,  e.g.,  the  stations  shall 
not  all  be  cclincar  and  N  =  2  (linear  case)  and  N  =  5  (quadratic  case). 


1 
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Appendix  A 

DERIVATION  OF  PARAMETER  ESTIMATE  EQUATIONS 

QUADRATIC  MODEL 

It  will  be  convenient  to  relabel  the  variables  of  Eq .  (1)  as 
fo 1  lows :  Let 


z(1)  -X.,  Z.(2)  -Y„  Z<3)  -X.  V.,  z<4>  =x.2,  Z.(i)  =v.2 

1  11  11  111  11  1 


Then  Eq .  (1)  can  be  written  in  matrix  form  as 


'  =  Z  A 


(A- 1 ) 


T  =  =  h  '  1  (column  matrix),  N  2  5 

A  =  -A. j  =5x1  (column  matrix),  of  unknown  parameters  A.  (A-2) 
Z  =  z(2\  Z(5);  =  N  v  5 

and  is  an  N  x  1  column  matrix,  u  =  1,  2,  ...,  5. 

•V 

In  particular,  values  of  A,  noted  as  A  ,  are  sought  which  minimize 
the  quadratic  form 


Q  =  (Z  A  -  T)/p(Z  A  -  T) 


(A-3) 


where  T  is  the  N  x  1  column  matrix  of  observations  of  T.  and  the 

l 

prime  indicates  the  transpose.  Upon  setting  the  gradient  Q  =  0  one 

*  (  1) 

obtains  the  well  known  result 


*  ,  - 1  / 

A  =  (Z  p  Z)  Z  p  T 


(A-4) 


* 


(  )  1  indicates  the  inverse  of  the  matrix  (  ),  and  1  the  transpose 


12 


where  it  is  assumed  that  the  columns  of  Z  are  linearly  independent  so 
that  (z'  p  Z)"1  exists. 

it 

It  is  easily  demonstrated  that  A  is  unbiased;  that  is 

E  (A*)  =  A 


(A-5) 


The  covariance  matrix  of  A  is  given  by  (see  Ea .  (A)) 


B*  =  e[[A*  -  A] [A*  -  A]'  I  =  /(Z'  P  Z)'1  Z7  p  V  p  Z(z'  p  Z)_1  Case  B 

0  /  -  -1  „  /  2  —  t  —  t  r.  \  1  - - A 


(Z7  p  Z)  z'  p  Z(Z7  p  Z )"  Case  A 


If  p  =  I,  the  matched  least- square  case  gives  (y  -  I) 


B*  »  a2(Z7  Z)"1 


(A-6) 


(A- 7) 


It  is  well  known  that  case  7d,  the  minimum  variance  estimator, 

(1) 


is  given  by 


A  =  (z'  r1  Z)"1  Z‘  if”1  T 


(A-8) 


and  the  corresponding  smallest  covariance  matrix  for  the  matched 
correlated  case,  corresponding  to  pY  =  I,  is 

S  «  a2(Z;  i'"1  Z)_1 


(A-9) 


LINEAR  MODEL 


The  derivation  for  the  linear  case  is  the  same  as  the  quadratic 
case  except  that  since  A3  =  A^  =  A^  =  0,  the  definition  of  Z  in 


z  -  [Z(1\  z(2)] 


Eq.  (A-2)  is  changed  to 


1 
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and  N  2  2  is  required, 
the  above  changes.  For 
of  5  x  5  and  A  is  a  2  x 


All  equations  (A- 1)  to  (A-9)  then  hold  with 
* 

example,  B  and  B  are  2x2  matrices  instead 


1  ins tead  of  a  5  *  1 • 
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Appendix  B 

DERIVATION  OF  CURVE  FITTING  EQUATION'S 

It  is  assumed  that  a  plane  acoustic  wave  is  incident  at  the 
array  at  bearing  angle  0  and  elevation  angle  ^  as  defined  in 

Fig.  1.  Since  the  quadrant  is  assumed  known,  there  is  no  loss  in 
generality  by  assuming  the  wave  as  incident  in  the  first  quadrant 
such  that 


o  *  t  -  0  •  e  •  f 

The  equation  of  the  phase  plane  is 

(sin  t>  cos  0)x  +  (sin  0  sin  0)y  +  (cos  b)z  -  P  =  0  (B-l) 

Consider  the  position  of  the  phase  plane  when  the  plane  is  incident 

at  the  reference  station  with  coordinates  (x  ,  v  ,  0) :  then  P  is 

o  o 

given  by 

P  =  sin  0  (cos  0)x  +  (sin  0)y 

o  o 

a nd  the  equation  of  the  phase  plane  is 

sin  t1  -(cos  0)  (x  -  x  )  +  (sin  0)  (y  -  y  )]  +  (cos  ())  z  =  0  (B-2) 

o  o 

It  is  required  to  compute  the  transit  time  of  the  phase  plane 

from  its  position  when  incident  at  station  (x  ,  y  )  to  the  time  when 

o  o 

the  phase  plane  is  incident  at  (x,  y)  .  Note  first  that  the 
distance  of  the  point  (x,  y)  from  the  phase  plane  through  station 
(xq,  y^ )  as  given  by  Eq .  (B-2)  is 


yr_ 
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1 1 - 
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P  =  -sin  0  (cos  0 )  (x  -  xq)  +  (sin  0)  (y  -  y  )]  (fi-3) 

where  xq  and  yQ  are  selected  such  that  -  (x.  -  x  )  ^  0,  -  (y .  -  y  )  s  0 

l  o  i  } o' 

for  each  of  the  station  coordinates.  The  transit  time  is  given  by 

(2  2) 

ray  theory  as  ’  ' 


where 


T 


n  (r) 


dr 


n(r) 


c  +  v  •  n 


c  +  v 


(B-4) 


(B-5) 


is  the  index  of  refraction  at  a  distance  r  along  the  ray  from  the 

station  at  (x,  y)  to  the  plane,  given  by  Eq.  (B-2),  formed  by  a  line 

perpendicular  to  the  plane  and 

Co  =  nominal  velocity  of  sound  =  344  m/sec  at  20°C 
— * 

n  -  unit  vector  in  the  direction  of  wave  propagation,  or 
perpendicular  to  the  phase  plane 

—♦ 

v(r)  =  wind  velocity  vector 
c  =  local  velocity  of  sound 
Vn  =  v  '  n  projection  of  v  on  n 

It  is  assumed  that  the  medium  is  horizontally  stratified  so  that 
both  c  and  v  are  functions  of  height  only.  For  a  standard  atmosphere 
one  may  write 


where 


c  (z) 


z 


Os  z  s  10  km 


dc 

dz 


4.4  meters/ sec/km 


(2) 


(B-6) 
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v(z)  versus  z  increases  logarithmically  with  z  for  heights  up  to  30 

(3) 

to  50  meters  and  then  at  a  slower  rate.  However,  for  the  purpose 

of  this  discussion,  wind  height  will  be  considered  a  slowly  increasing 
linear  function  of  height  which  can  be  represented  over  the  range  of 
altitudes  of  interest  as 


v  (z)  =  v  (o)  +  K.  z 
n  n  J. 

Setting  z  =  r  cos  0,  Eq .  (B-4)  becomes 


(B-7) 


v  (z) 


c  cos  0  L  c  dz  c  J 


dz 


(B-8) 


where  p  is  sufficiently  small  so  that  ;  ~[^/c0  |  <  <  and 

(v  ( z ) / c  )  <  <  1.  Substituting  Eq.  (B-7)  into  Eq.  (B-8)  and  integra- 
n  o 

ting  Eq.  (B-8)  gives 


r  v  (° ) . 

t  « M  i  -  — — 

c  L  C_ 


1  -  i2-) 

cos  t> 

+  Ki]" 

\c  / 

J  o 

2 

Idz  1J 

Equation  (B-9)  is  a  quadratic  in  p  which  can  be  written  in  the  form 

.  2 
t  »  a  P  +  p  p 

On  substituting  Eq .  (B-3)  into  Eq .  (B-9)  one  finds 

2  2 

t  =  A,  X  +  A„  Y  +  A,  X  Y  +  A,  X  +  A.  Y 


(B-9) 


(B-10) 


(B- 11) 


The  coefficients  are  given  by 
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-  (y  sin  fa  cos  0)  d 

A-,  =  (y  sin  fa  sin  0)  d 

2  9 

=  (2?  sin"  fa  sin  G  cos  0)  d“ 

2  9  0 

A^  =  (3  sin"  fa  cos"  3)  d" 

2  2  9 

A^  =  (3  sin"  fa  sin"  0)  d" 


O' 


L 


1  - 


vn<°>. 


c 

o 


3 


cos  fa  dc 

0  Id/.  4 


K 


iJ 


Define  the  effective  ground  trace  velocity 


c  by 
8 


c  =  ot  sin  fa 
8 


(B- 12) 


(B- 13) 


Then  estimates  of  both  c  and  9  may  be  obtained  as  follows: 
Note  that 


(B— 14) 

(B-15) 


Thus  the  estimates  of  A^  and  A0  provide  estimates  of  the  bearing 

angle  0  and  the  effective  ground  trace  velocity.  The  quantity  t  is 

o 

the  time  for  the  wave  to  travel  a  distance  d  on  the  ground. 

If  3  =  0,  so  that  the  linear  model  for  T  can  be  used,  the  amount 
of  data  processing  is  reduced  and  the  rms  of  the  estimates  of  0  and 
Tq  is  decreased.  From  Eq.  ( B -  1 1 ) 
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0/ 


p  V  + 1 P1 


From  Eq .  (B-12) 


2  2 

A1  +A2 


d  sin  0 


A4  +  A5 

2  2 
d“  sin"  0 


so  that  Eq.  (B- 16)  can  be  written 


=  O'  P 


1  + 


A.  +  A 
4 _ 

2  2~ 

A1  +A2 


d  sin  0 


The  value  of  p/ d  sin  0  is  clearly  determined  from  Eq .  (B-3)  as 


p/d  sin  t  =  X  cos  »  +  Y  sin  9 

so  that  the  maximum  magnitude  of  p/ d  sin  t  =  the  maximum  normalized 
dimension  of  the  array.  Let  this  characteristic  value  be  D  where 


D  =  max 


d  sin  9- 


0,  V 


Then ,  i 1 


A4  +  A5 


D  <  <  1 

v. 

2 


one  may  take  0  =  0,  and  therefore  A^  -  A^  =  A^  =  0,  and  use  the 


(B  - 1 6) 


(B-17) 


(B- 18) 


linear  model . 


20 


Appendix  C 

BEARING  ANGLE  ACCURACY 

The  relationship  between  the  bearing  angle  0  and  the  coefficients 
of  the  curve  fit  A^  ,  A0  is  given  in  Eq .  (B-14).  This  relationship  is 
nonlinear.  However,  if  the  errors  in  the  coefficients  are  small, 
then  the  errors  in  0  can  be  determined  as  follows 


tan  0  =  A^/A 


A0  =  cos  © 


1  {A2fiAl  -  A  4A,} 


(C-1) 


where  £0  is  the  random  error  in  0  due  to  random  errors  AA^  and  £A^  in 
the  parameter  estimates  A^  and  \^.  Since  E(At\)  =  0,  then  E(AA.)  =  0, 
j  =  1,  2,  3,  ,  5,  and  E(A0)  =  0.  Then 


,  4 

cos  0\ 


E(A0  )  =  [A,  -  Ax]  B2  [A2  -  Aj]' 


where  B2  is  the  2X2  submatrix  of  the  covariance  matrix  given  by 
Eq.  (A-7)  or  Eq.  (A-9);  e.  g. 


B11 

B12 

2 

bn 

b12 

B„„ 

=i  c 

h  , 

b 

21 

22 

21 

22 

B2  = 


where  the  b. .  are  the  normalized  covariance  E(AA.  AA . ) 

ij  i  J 


(C-2) 


(C-3) 

(C-4) 


b.  .  SI 
ij  1 


The  value  of  A^  used  in  the  estimation  is  matched  to  the  appropriate 
choice  of  p  for  a  given  i|r  to  determine  which  B  matrix  to  use. 

From  Eq.  (B-12) 
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sc  that 


A  =  T  cos  9,  A  =  T  sin  9, 
1  o  2  o 


E(A9  )  = 


f)  { 


cos  9  b, ,  +  sin  9  b„,,  -  2  b,.  cos  9  sin  9 


11 


22 


12 


(C-5) 


(C-6) 


Define 

L2(9)  =  (C-7) 

k'V 

the  normalized  variance  of  9. 

2 

Equation  (C-6)  shows  that  E(A9  )  is  inversely  proportional  to 
T  ,  the  ground  transit  time  of  the  wave  over  the  distance  given  by 
the  scale  factor  d.  Assuming  d  to  be  fixed  (say  the  x,  y  coordinate 
distance  between  adjacent  stations  in  a  square  array),  then  -*  0 
as  0  -*  0.  (See  Fig.  1.)  In  this  case  the  ground  trace  velocity  is 
infinite  and  9  becomes  indeterminant,  as  is  expected.  Thus,  it  is 
required  to  limit  0  so  that  0  2:  0q  before  an  attempt  to  estimate  9 
is  considered.  Define 


c.q  *\/E(A92)  =  ?-  ’  L(9)  in  radians 

o 

L(9)  is  shown  in  Tables  1  and  2  and  gives  the  bearing  angle  accuracy 
in  radians . 

From  Eq.  (C-6)  note  that  if  b^  =  b ^ 

E(A92)  =  (^—J  [b^  -  2  b12  cos  9  sin  9}  (C-8) 

o 

so  that  L(9)  is  symmetrical  with  respect  to  9  *  45°.  When  the 


stations  are  placed  symmetrically  with  respect  to  the  line  y  =  x, 


NORMALIZED  BEARING  ANGLE  rms  (L(^))  AND  RATIO3  R  VERSUS 
CORRELATION  PARAMETER  k  FOR  THE  LINEAR  CASE 


y 

i 


i 
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f  1  is  the  number  of  stations 
is  the  bearing  angle 


NORMALIZED  BEARING  ANGLE  rms  (L(0))  AND  RATIO  R  VERSUS 
CORRELATION  PARAMETER  k  FOR  THE  QUADRATIC  CASE 
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then  it  is  obvious  that  one  may  interchange  y  and  x  and  demonstrate 
that  bu  =  b22  so  that  the  symmetry  conditions  given  by  Eq .  (C-8) 
hold . 

Finally,  Eq  •  (C-l)  seems  to  require  0  •*  rr/2.  However,  one  can 
define  cotan  =  A^/A2  and  derive  Eq.  (C-6)  as  the  end  result  so  that 
Eq .  (C-6)  does  hold  for  all  0. 

TRANSIT  TIME  ERROR 
From  Eq.  (B-15) 

..  ~A1  +A2  AA2 

A  —  ———————— 

o  T 

o 

so  that 

E(At  )  =  0 
o 

2  2'  2  2 

E(Ato  )  =  -  '.b^  cos  ®  +  ^22  s^n  ®  +  2cos  6  sin  b^2i 

If  bll  =  b22’  wbicb  bolds  for  stations  symmetrically  placed  with 
respect  to  the  line  Y  =  X 

E(AT  2)  r  2 

- — —  =  [b^  +  2  cos  0  sin  0  b^}  =  M  (0) 

rr 

The  normalized  variance  M^(0)  is  symmetric  with  respect  to  0  =  45°. 
Note  that  when 

0  =  0  or  0  =  tt/2,  M(0)  =  L(0) 

o 

for  any  value  of  0 

M  (9)  -  L  (0)  =  4b^2  sin  9  cos  0 


(C-9) 


(C-10) 


(C-ll) 


(C-12) 


(C-13) 


and  therefore 
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M2^)  -  L2(9)  S  2  b12 

the  equality  sign  holding  for  9  =  45°  when  the  symmetry  conditions 
b^  =  b^2  hold  and 

M2  (9)  +  L2(9)  =  2  bn 

independent  of  0. 

Table  3  presents  M(9) ,  9  =  45°,  for  the  linear  and  quadratic 
case  with  k  from  .125  to  256  and  values  of  N  as  indicated. 


(C- 14) 


(C- 15) 


i 
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Appendix  D 

COMPUTATION  OF  NORMALIZED  BEARING  ACCURACIES 
L(°)  and  M(9)  FOR  A  SQUARE  ARRAY 

In  this  appendix  computations  of  L(P) ,  the  normalized  bearing 
accuracy,  and  M(P),  the  normalized  ground  trace  timing  accuracies  given 
by  Eq.  (C- 1 1 )  ,  are  described  for  a  specific  array  configuration  shown 
in  Fig.  2. 

The  numbers  in  Fig.  2  show  the  normalized  coordinates  (X,  Y)  and 
station  index  number. 

For  the  linear  case,  L(Q)  and  M(9)  are  computed  for  the  first 
four  stations  (N  =  3),  the  first  eight  stations  (N  =  7),  which  includes 
the  previous  stations,  and  all  the  stations  (N  =  15). 

For  the  quadratic  case,  L(P)  is  computed  for  N  =  7  and  N  =  15 
defined  over  the  same  set  of  stations  as  in  the  linear  case  for 
corresponding  N. 

Computations  are  performed  for  9=0°,  15°,  30  and  45  for 
values  of  the  correlation  parameter  k  =  (2)^,  j  =  -3  to  8,  in  steps 
of  1.  For  small  values  of  k  (.125  or  .25)  the  effect  is  essentially 
the  same  as  taking  ilr  =  I,  so  that  this  case  will  not  be  computed  sepa¬ 
rately.  For  large  values  of  k,  the  At  errors  at  each  station  are 
heavily  correlated,  and  one  may  note  the  effect  of  using  a  mismatched 
processing  such  as  least  squares  on  this  data  versus  using  the  matched 
processing,  p  v  =  I,  of  Eq.  (7d).  The  Fortran  program  from  which  L(0) 
and  M(9)  are  computed  is  shown  on  the  following  pages.  Figures  3  to 
12  present  the  L(9)  values  graphically  for  possible  interpolation  and 
visual  comparison.  Tables  1  to  3  present  numerical  results  of  the  program. 


Note  that  the  Fortran  program  is  sufficiently  general  to  handle 


an  arbitrary  set 
above,  provided  N 
case,  and  not  ail 


of  stations  and  not  just  the  square 
^  2  in  the  linear  case  and  N  i  5  in 
the  stations  are  colinear. 


array 

the 


described 

quadratic 
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C 

C  PROGRAM  TO  COMPUTE  NORMALIZED  BEARING  ACCURACY  AND  NORMALIZED  GROUND 
C  TRACE  TIMING  ACCURACY  FOR  LEAST  SQUARES  AMD  'MNIMIJM  VARIANCE  CASES. 
REAL  K 

DIMENSION  ZlM 5)»Z2<35).Z3<35),Z4(35)  ,75(35)* NN (35)  ♦  C  A  Y ( 3  5  )  ♦  T  MET ( 3 
15 ) *S( 1 *1 ) *Z ( 35  *35)  *R( 35*35) *SINE ( 35 ) ,C0SINF( 35 ) ,PSI ( 35*35)  *P( 35*35 
2)*ZT(35*35)*ZTZ(35, 35).  Zll ( 35  *35 )  *  712(3 5*35)*Z13(35*35)*B(35*35)*Z 
321(35*35)  *C( 35*35) *CT( 35*35) *Z31 ( 35*35) »EL£Q(35*35) *EL ( 3  5  )  *Z 41(35. 
43  5 ) *EMSQ( 35 • 35 )  *FM( 35 )  *  IP  I VOT ( 35 )  .  1 N^EX (^5*2) *ELL ( 35 ) *  EMM ( 35)  *  F  L  R ( 
535) » EMR ( 3  5 ) 

READ  2*N2*LB2.LC?*L  12*LA2*LF2 
READ  1  » ( Z 1 ( N )  *N  =  1 *N2  ) 

READ  1  » ( Z2 (N )  *N*1 *N2  ) 

READ  1  *  <  Z  3 ( N )  *N  =  1 *N2) 

READ  1*(Z4(N) *N*1 *N2) 

READ  1  * ( Z  5 ( N )  »  N  =  1  ♦  N  2 ) 

READ  2  ♦  ( NN ( LB ) ,LB=1 *LB2 ) 

READ  3.  ( C A Y ( LC )  *LC  =  1  *LC2 ) 

READ  4  ,  (  T  H  E  T (LI  )  *  L  I  = 1  *  L I  2 ) 

1  FORMAT ( 18F4.0 ) 

2  FORMA  T ( 1 8  I  4 ) 

3  FORMAT  I 8F9.3 ) 

4  FORMAT ( 12F6* 2 ) 

S  (  1  ,  1  )  *  1  . 

Z  MATRIX  (N?  X  5)  IS  EORvrn. 

DO  10  1=1,5 
DO  10  N= 1 ,N2 

IF  (  I  . EO . 1 )  Z (N  ♦  I  )*Z1  (  N) 

IF  (  I.E0.2)  Z  (  N  *  I  )=Z2(M 
IF  (  I  .EQ.3)  Z  (N.I  )  —  Z 3 ( N ) 

IF  (  I.FQ.4)  Z ( N  »  I  )=Z4(N) 

10  IF  (  I  .  EQ . 5 )  Z ( N  ,  I  ) =25 ( N) 

C  R  MATRIX  (N?  X  N2 )  IS  FORMED. 

DO  20  1=1, N2 
DO  20  N= 1 »  N2 

20  R (  I  » N ) =  SQRT(  ( Z  1  (  I  )  -2 1 ( N )  ) **2  +  <  Z  2  (  I  )-72(M)  ) **2 ) 

C 

C  SINE  AND  COSINE  VALUES  ARE  CALCULATED  HERF  TO  SAVE  TIME. 

RAD=1 .74532925E-2 
DO  25  L I  *  1 , L I  2 
THETA*THET ( L I )  *RAD 
S I NE ( L I ) *S I N  (  THETA) 

25  COSINE(LI ) »COS( THETA) 


n  n  no  noon 
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PROBLEM  BEGINS. 

LINEAR  WHEN  LA=1.  QUADRATIC  WHEN  LA  =  2. 

DO  100  L A  = 1 »  L  A2 
IF  (LA.NF. 1 )  GO  TO  27 
LB  1  =  1 
M  =2 

GO  TO  28 

27  LR 1 =2 
M  =  5 

A  VALUE  OF  N  IS  PICKED. 

28  DO  lOu  LB=LB1.LB2 
N  =  NN( LB  ) 

A  VALUE  OF  K  IS  PICKED,  PSI  MATRIX  (N  X  N )  IS  FORMED. 

DO  100  LC  = 1 ♦ LC2 
K  =  C AY ( LC  ) 

IF  (LA.FO.l)  PRINT  2000,N,K 
IF  (LA.EO.2I  PRINT  2001, N»K 
2000  F0RMAT(1H1,2X,6HLINEAR,4X,2HN«I2,4X,2HK*F8.3///) 

200  1  FORMAT ( 1 H 1  , 2X.9HQUADRATIC ,4X, 2HN=  12 »4X , 2HK»FB . 3/// ) 

DO  30  LE  =  1  ,N 
DO  30  LD= 1  ,N 

PSKLD.LE)  =EXP(-R(LD,LE1/K) 

30  P(LD,LE)-PSI (LD.LE) 

C  CASE  1.  RHO  MATRIX  =  IDENTITY  (N  X  N)  IS  FORMED.  (MISMATCHED) 

C  ZT  MATRIX  (M  X  N)  =  Z  (N  X  M)  TRANSPOSE  IS  FORMED. 

C  ZTZ  MATRIX  (M  X  M)  =  MATRIX  PRODUCT  OF  ZT  AND  Z  15  FORMED. 

C  ZTZ  INVERSE  MATRIX  (M  X  M)  IS  FORMED. 

C  B  (NORMALIZED  COVARIANCE  MATRIX  -  M  X  M)  *  MATRIX  PRODUCTS  OF 
C  ZTZ  INVERSE  (M  X  M),  ZT  (M  X  N),  PSI  (N  X  N),  Z  (N  X  M),  ZTZ  INVERSE 
C  (M  X  M)  IS  FORMED. 

DO  80  LF* 1 , LF 2 
IF  (LF.NE.  1  )  GO  TO  50 
DO  40  I  *■*  1  »  M 
DO  40  J* 1  ,  N 
40  ZT ( I  *  J ) *Z ( J  •  I  ) 

CALL  MATMUL  ( ZT , M , N , Z , M , Z T Z  ) 

CALL  MAT  I N V  ( ZT Z ,M , S ,0 , 1 P I VOT , I NDEX , I S I NG ) 

IF  ( ISING.NE.O)  GO  TO  71 
CALL  MATMUL  ( ZTZ , M , M ,ZT , N » Z 1 1  ) 

CALL  MATMUL  ( Z 1 1  ,M ♦ N , PS  I , N , Z1 2 ) 
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CALL  MATMUL  ( Z12.M.N.Z.M.Z13) 

CALL  MATMUL  (Z13.M.M.ZTZ.M.B) 

PRINT  1002 

1002  FORMAT (  1H0.4X .1 RHCASE  1  ( M I SMATCHFD  )  ) 

PRINT  1003 

1003  FORMAT ( //7X * 32HR  (NORMALIZED  COVARIANCE  MATRIX)/) 

DO  40  1  =  1  .M 

43  PRINT  1004  * ( B ( I  * J ) * J* 1 *M  ) 

1004  FORMAT ( 1H  .  F 1 7 . 8 *4F 15. 8  ) 

GO  TO  60 

C 

C  CASE  2.  RHO  MATRIX  (N  X  N)  =  PSI  INVERSE.  (MATCHED) 

C  PSI  INVERSE  MATRIX  (N  X  N)  IS  FORMED. 

C  B  (NORMALIZED  COVARIANCE  MATRIX  -  M  X  M)  a  MATRIX  PRODUCTS  OF  ZT 
C  N ) *  PSI  INVFRSE  (N  X  N).  Z  (N  X  M)  INVERSE  IS  FORMED. 

50  CALL  MAT  I NV  ( P * N * S *  0  *  I P I VOT ♦ I NDEX ♦ I S I NG ) 

IF  ( ISING.NE.O)  GO  TO  72 

CALL  MATMUL  ( Z T *M * N *P ♦ N * Z 2 1  ) 

CALL  MATMUL  ( Z ? 1 ♦ M ♦ N . Z . M ♦ B ) 

CALL  MAT  I N V  ( B .M * S .0 » I P I VOT  ♦  I NDEX » I S I NG ) 

IF  ( ISING.NE.O)  GO  TO  73 
C 

C  NORMALIZED  COVARIANCE  MATRIX  B  IS  PRINTfD  OUT. 

PRINT  1005 

1005  FORMAT ( //AX » 16HCASE  2  (MATCHED ) > 

PRINT  1003 

DO  51  1  =  1  .M 

51  PRINT  100A. ( B( I »J) »J*1 »M) 

C 

C  FOR  MISMATCHED  AND  MATCHFD  CASES.  NORMALIZED  BEARING  ACCURACY  FL 
C  NORMALIZED  GROUND  TRACE  ACCURACY  EM.  AND  PATIOS  OF 
C  LR  *  EL  (MISMATCHED)  /  EL  (MATCHFD)  AND 
C  MR.*  EM  (MISMATCHED)  /  EM  (MATCHED)  ARE  PRINTED  OUT. 

60  DO  62  L I  =  1  .  L  I  2 
C(  1  .1 ) “COSINE (L  I  ) 

C(2.1)=-SINE(LI  ) 

C( 3.1 ) *0. 

C  (  4 . 1 ) *0 . 

C( 5*1 ) *0. 

CT ( 1 . 1 ) =COS I NE ( L I ) 

CT ( 1 »  2 ) *-S I NE ( L I ) 

CT  (  1  *  3  )  =0 . 

CT ( 1 .4 )*0. 

CT (  1 .5 )*0. 

CALL  MATMUL  ( CT . 1 *M * B . M » Z 3 1 ) 


(M  X 


33 


CALL  MATMUL  ( Z 3 1 » 1  * M ♦ C . 1 . EL SQ ) 

FL ( L I )=SQRT(CLSQ( 1.1) ) 

C  (  2  ♦  1  )  =  S  I  N  F  (  L  I  ) 

CT (1 ,2 )  =  S I  ME  < LI  ) 

CALL  MATMUL  (  CT .1 »M.B.M.Z41  ) 

CALL  MATMUL  ( Z4 1 . 1 ♦ M ♦ C . 1 . EMSQ ) 

F M ( L I  )=SQRT(EMGQ( 1.1  )  ) 

IF  (LF.NE.l )  GO  TO  62 
FLL ( L I )  =  EL (LI ) 

EMM ( L I  )  =  E  M ( L I  ) 

62  PRINT  1006. THET  <  L I ) ♦  FL (LI ) ♦  THET (LI ) »FM(LI  ) 

1006  FORMAT  (/(7X.?HL(F6.?.2H)=F1?.8.4X.?HV(F6.2.?H)=F12.8) ) 
GO  TO  80 

71  PRINT  1007 

1007  FORMAT ( 1MC .4X  .24HZTZ  INVERSE  IS  SINGULAR.) 

GO  TO  80 

72  PRINT  1008 

1008  FORMAT ( 1HC.4X.24HPSI  INVERSE  IS  SINGULAR.) 

GO  TO  80 

73  PRINT  1009 

1009  FORMAT ( 1H0.4X.30HB  MATRIX  (CASE  2)  IS  SINGULAP.) 

80  CONTINUE 
PRINT  81 

81  FORMAT ( 1H0 ) 

00  90  L I  =  1  .  L  I  2 
FLR(LI )  =  ELL(LI  )  / E  L ( L  I  ) 

EMR(Ll  )  =  E  M  M ( L I ) / E  M  <  L I  ) 

90  PRINT  101C .THET (LI ) .ELR( L  I  )  .THET ( LI ) .EMR( LI  ) 

1010  FORMAT ( 1H  ,4X  »3HLR(  F6.2 .2H) *F 1 2 • 8 .4X.3HMR ( F6.2 .2H)  =  F  1 2 . 8  ) 
100  CONTINUE 

CALL  EXIT 
END 


i.  -4 
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C  MATRIX  MULTIPLICATION 

SUBROUTINE  MATMUL  ( A ♦ M » N  ♦  R » L ♦ C  ) 
DIMENSION  A(35.35)»R(35*35)*C(35.36) 
DO  20  I ■ 1  * M 
DO  20  1 ♦ L 

SUM=0. 

DO  10  J *  1  *  N 

10  SUM*SUM+a(  I  *  J ) *  B ( J*<) 

20  C  (  I  ♦  K  )  «SUM 
RETURN 
END 


C  MATRIX  INVERSION  WITH  ACCOMPANYING  SOLUTION  OF  LINEAR  EQUATIONS 

SUBROUTINE  MAT  I  NV  (  A  .  N  .  B  .  M  .  I  P I  VOT  ♦  INDEX  ♦  IS  INC.) 

C 

DIMENSION  A<35*35)*B(1*1) ♦  I  PI VOT ( 35) . INDEX! 35*2) 

EQUIVALENCE  (IR0W*JR0W)«  l I COLUM * JCOL UM )  ♦  (T*  SWAP) 

C 

C  INITIALIZATION 

5  ISING  *  0 
15  DO  20  J* 1 ♦ N 
20  I P I VOT ( J ) *0 
C  BIG  LOOP  ON  I 

30  DO  550  I  a  1  * N 
35  IROW  »  0 
40  AMAX-0.0 

C  SEARCH  FOR  PIVOT  ELEMENT 

45  DO  105  J* 1 ♦ N 

50  IE  (  I P I VOT ( J ) # EQ • 1  )  GO  TO  105 

60  DO  100  KM,N 

70  IE  (  IPIVOT ( K ) *EO* 1  )  GO  TO  100 

BO  IF  (ABS(AMAX).GE.ABS(A(J»K) )  )  GO  TO  100 

05  IROW»J 

90  ICOLUM-K 

95  AMAX«A(J,K) 

100  CONTINUE 
105  CONTINUE 

107  IF  ( I  ROW*  EQ*  0 )  GO  TO  750 


•* 


V 


n  n  r>  ^  n  n  nnn  n  n  r> 
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110  I P I  VO  T (  ICCIUM)  =  1 

INTERCHANGE  ROWS  TC  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

130  II-  (  IROW.EQ.  ICOLIJM  )  GO  TO  260 
150  DO  200  L  =  1  » N 
16  '  S  W  A  P  =  A  (IROW.L) 

170  A (  I  ROW  » L ) =  A (  I  COL  UM  »  L ) 

200  A ( I COLUM.L ) *SWAP 
205  IF  (M.LF.O)  GO  TO  260 
210  DO  260  L  =  1  *  M 
220  SWAP=B( IROW.L) 

230  B( IROW.L )«B( I COLUM.L) 

250  B ( ICOLUM.L )=SWAP 
260  INDEX ( I ,1 ) = IROW 
27^  INDFX ( I ,2 ) = ICOLUM 

OIVIDF  PIVOT  ROW  HY  PIVOT  FLFMENT 

33^  A(  ICOLUM,  irOLlJM)=1.0 
340  DC  360  L  =  1  *  N 

350  A ( ICOLUV,L)=A( I COLUM.L ) /AMAX 
355  IF  (M.LF.O)  GO  TO  380 
360  DO  370  L*1 .m 

3  70  B  (  I  COLUM.L  )  =13  (  I COLUM. L ) / AMAX 
COMPLETF  THE  PIVOT 
380  DO  550  L1*1*N 

390  IF  (LI .FO. ICOLUM)  GO  TO  550 
400  T=A (LI . ICOLUM) 

420  A ( LI . ICOLUM) =0.0 
430  DO  450  L  =  1  » N 

450  A ( L 1 ;  L ) = A ( L 1  *  L ) “A ( ICOLUM,L)*T 
455  It  (M.LF.O)  GO  TO  550 
460  DO  500  L= 1 *M 

500  B(L1.L)=B(L1.L)-B( ICOLUM.L)*T 
550  CONTINUE 

INTERCHANGE  COLUMNS 

600  DO  710  I«1 «N 
610  L  =  N+ 1 - I 

620  IF  (  INDEX(L.1  )  . EO. I NDE X ( L . 2 )  )  GO  TO  710 
63C  JROW= I NDE X ( L  » 1  ) 
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640  JC0LUM=INDFX(L.2) 

650  DO  705  K* 1 *N 
660  SWAP*A ( K » JROW ) 

670  A ( K » JROW) «A ( K . JCOLUM ) 
700  A(X»JC0LUM)=SWAP 
705  CONTINUF 
710  CONTINUF 
740  RFTURN 
C  SINGULARITY  FLAG 

750  ISING  =  1  +  N  -  I 
760  GO  TO  740 
END 


* 
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INPUT 


* 


Input 

Order 

Variab  le 

Name 

Explanation  of  Variable 

Restriction 

Input 

Format 

1 

N2 

Total  number  of  stations 

£ 

36 

1814 

1 

LB  2 

Number  of  sets  of  stations 

1814 

1 

LI2 

Number  of  0  values 

35 

1814 

1 

LC2 

Number  of  k  values 

£ 

35 

1814 

1 

LA2 

Linear  case  only  when  LA2=1, 
linear  and  quadratic  case 
when  LA2=2 

1814 

1 

LF2 

Mismatched  case  only  when  LF2=1, 
mismatched  and  matched  cases 
when  LF2=2 

1814 

2 

Z1 

x  coordinate  of  each  station 

£ 

35 

18F4.0 

3 

Z2 

y  coordinate  of  each  station 

£ 

35 

18F4.0 

4 

Z3 

Product  of  x  and  y  coordinates 
of  each  station 

18F4.0 

5 

Z4 

x2  of  each  station 

18F4.0 

6 

Z5 

2  , 

y  of  each  station 

18F4.0 

7 

NN 

Number  of  stations  minus  one 
used  in  each  set 

NN>2  for 

1  inear 

NN">5  for 
quadratic 

1814 

8 

CAY 

Correlation  parameter 

0<k< 

10,000 

8F9.3 

9 

THET 

Bearing  angle  of  plane  wave 
in  degrees 

12F6.2 

The  program  is  designed  to  handle  the  parameters  as  indicated  by  the 
restrictions.  The  program  however  has  not  been  checked  to  the  limit 
of  these  restrictions.  Computations  for  the  cases  presented  in  this 
report  have  been  checked.  Other  cases  may  require  additional  verifica¬ 
tion. 
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Sample  Input 


15  3  12  4  2  2 

3.  3.  0.  1.  3.  2.  3. 

0.  3.  3.  3.  1 .  3.  2. 

0.  9 .  0 •  3 •  3 •  6 •  5* 

9.  9.  0.  1.  9.  4.  9. 

0.  9.  9.  9.  1.  9.  4. 

3  7  15 

.12-5  .25  .5 

32.  64.  128. 

0.  15.  30.  45. 


2.  1. 
2.  2. 
4.  2. 

4  .  1  . 

4  .  4  . 

1  . 

256. 


N2 

2.  1.  0.  2.  1.  0. 

1  .  1.  2*  0.  0.  1  • 

2.  1.  0.  o.  0.  0. 

4.  1.  0.  4.  1.  0. 

1  .  1.  4.  0.  0.  1* 

2.  4.  6. 


LB 2  LC  2  L I  2  L  A 2  LF? 
Z  1 
Z  2 
Z  3 
Z  4 
Z  5 
N 

16.  CAY 
CAY 
T  H  F  T  A 


OUTPUT 


N  =  Number  of  stations  minus  one 
K  -  Correlation  parameter 
L(6)  =  Normalized  bearing  angle  accuracy 


M(6) 


LR(6) 


Normalized  ground  trace  timing  accuracy 
L(b) 

_ misma  tched 

L(c)  ,  , 

matched 


MR(  6) 


M(3) 

W) 


mismatched 

matched 
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Sample  Output 


UNCAP  \l=  3  <=  1.000 


CA cr  1  (MISMATCHED) 

B  (NORMALIZED  COVARIANCE  MATRIX  ) 


0. 

07459377 

-0.03492071 

-0. 

03492071 

0.07459377 

L  ( 

o 

. 

ii 

0.27311860 

M< 

0  .  )  = 

0.27311860 

L  ( 

1  5 . 0  0  )  = 

0.30340423 

M  ( 

15.00)= 

0.23902597 

L  ( 

30.  l.O  )  = 

0.32378387 

M  ( 

30.00 ) = 

0.21059807 

L  ( 

45 .0C  )  = 

0.33092972 

M  ( 

45. PA )  = 

0. 19918097 

CASE  2  (MATCHED) 

P  (NORMALIZED  COVARIANCE  MATRIX) 
0.07457582  -0.03493867 


-0 

•  03493867 

0.07457582 

L  ( 

0 

.  )  = 

0. 27308570 

M  ( 

0.  )  = 

0.27308573 

L  ( 

15 

.00  )  = 

0.30338944 

M  ( 

15.00  )cF 

0.23896962 

L  ( 

30 

•  00  )  = 

0.32378016 

M  ( 

30.00 ) = 

0.21051850 

L  ( 

45 

.GO  )  = 

0.33092973 

M  ( 

45.00)* 

0.1 9909079 

LR  ( 

0. 

)  = 

1.00012037 

MR  ( 

II 

• 

o 

1.00012037 

LR  ( 

15. 

00  )  = 

1 .00004874 

MR  ( 

15.00 ) = 

1.00023580 

LR( 

30. 

00  )  = 

1 .00001143 

MR  ( 

30.00 ) = 

1.00037797 

LR( 

4  5 « 

00  )  = 

0.99999996 

MR  ( 

45.00 ) = 

1.00045294 

L  (9),  normalized  rms  bearing  angle  accuracy 
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9,  bearing  angle  in  deg 

Fig. 3— Normalized  rms  bearing  angle  accuracy  versus  bearing  angle 


1(0)/  normalized  rms  bearing  angle  accuracy 
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L(0),  normalized  rms  bearing  angle  accuracy 
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i 

1 

J 

*! 

1 

_ L 

"  r—  t  — i - 1 - 

- 1 - 1  j - 1 - 

- 1 - 1 - 1 - 1 - 

Matched  case 

P'P  ~  I 
k  =  0.125 

- 

- quadratic 

case 

linear  case 

- 

N  +  1  =  number  of  _ 
stations 

15 _ _ 

■ - 

— 

N  -  3 

- 

_ 7 _ 

15 

_ 1 _ 1 _ 1 _ 1 _ 

_ 1 _ 1 _ 1 _ 1 _ 

i _ 1 _ 1 _ 1 - 

0  15  30  45 

6,  bearing  angle  in  deg 


Fig. 6— Normalized  rms  bearing  angle  accuracy  versus  bearing  angle 


L(0),  normalized  rms  bearing  angle  accuracy 
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Fig. 7— Normalized  rms  bearing  angle  accuracy  versus  bearing  angle 


L(0),  normalized  rms  bearing  angle  accuracy 


45 


0,  bearing  angle  in  deg 


Fig. 8—  Normalized  rms  bearing  angle  accuracy  versus  bearing  angle 


smotched  case 


Normalized  •'ms  bearing  angle  accuracy  versus  correlation  parameter 


Fig.  10— Normalized  rms  bearing  angle  accuracy  versus  correlation  parameter 


k,  correlation  parameter 

Fig.  11— Normalized  rms  bearing  angle  accuracy  versus  correlation  parameter 
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